home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / MATH / MFLOAT10.ZIP / BESSEL.CPP < prev    next >
C/C++ Source or Header  |  1993-04-28  |  2KB  |  88 lines
  1. /* *** This test program shows the advantages of mfloat numbers. *** */
  2. /* *** The bessel function J0(x) is calculated by the series.    *** */
  3.  
  4. extern unsigned _stklen = 40000U;
  5.  
  6. #include <conio.h>
  7. #include <iomanip.h>
  8.  
  9. #include <math.h>
  10. #include "mfloat.cxx"
  11.  
  12. /*-------------------------------------------------------------------------*/
  13.  
  14. long double J0(long double x) {
  15.  
  16.   const mfloat epsi = 1e-20;
  17.   mfloat sum, sqrx, prod, mepsi;
  18.   int    i;
  19.  
  20.   sqrx = x * x / 4;
  21.   sum = 1;
  22.   prod = 1;
  23.   i = 0;
  24.   do {
  25.     i++;
  26.     prod = prod * sqrx / i / i;
  27.     if (i & 1)
  28.       sum -= prod;
  29.     else
  30.       sum += prod;
  31.   } while (epsi <= prod);
  32.   return mftold(sum);
  33. };
  34.  
  35. /*-------------------------------------------------------------------------*/
  36.  
  37. long double J0x(long double x) {
  38.  
  39.   const long double epsi = 1e-20;
  40.   long double sum, sqrx, prod;
  41.   int    i;
  42.  
  43.   sqrx = x * x / 4;
  44.   sum = 1;
  45.   prod = 1;
  46.   i = 0;
  47.   do {
  48.     i++;
  49.     prod = prod * sqrx / i / i;
  50.     if (i & 1)
  51.       sum -= prod;
  52.     else
  53.       sum += prod;
  54.   } while (epsi <= prod);
  55.   return sum;
  56. };
  57.  
  58. /* ------------------------------------------------------------------------ */
  59.  
  60. int main() {
  61.  
  62. int accuracy;
  63. long double x;
  64.  
  65.   clrscr;
  66.   cout << "              Calculation of the bessel function J0(x)\n"
  67.        << "              ========================================\n"
  68.        << "\n";
  69.  
  70.   x = 100;
  71.   cout << "mantissa words     result     (x = "
  72.        << x << ")\n\n";
  73.   for (accuracy=1; accuracy<=15; accuracy++) {
  74.     setmantissawords(accuracy);
  75.     cout << setw(3) << accuracy
  76.      << "                J0(" << setw(5) << setprecision(1)
  77.      << x << ") = "
  78.      << setiosflags(ios::scientific)
  79.      << setw(25) << setprecision(18) << J0(x) << "\n"
  80.      << resetiosflags(ios::scientific);
  81.   }
  82.   cout << "\nIEEE arithmetic    J0(" << setw(5) << x << ") = "
  83.        << setiosflags(ios::scientific)
  84.        << setw(25) << setprecision(18) << J0x(x) << "\n"
  85.        << resetiosflags(ios::scientific);
  86.   cout << "\nThe accuracy of 12 mantissa words is needed to get a good result!\n";
  87.   return(0);
  88. }